Data analysis for:
Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.
Prokaryotes Become Larger at High Temperatures but They Do Not Grow Faster
Library and Session info
library(AICcmodavg)
library(ape)
library(geiger)
library(phytools)
library(phylolm)
library(lattice)
library(MuMIn)
library(nlme)
library(raster)
library(rphylopic)
library(scales)
library(vioplot)
sessionInfo()> R version 4.4.1 (2024-06-14)
> Platform: aarch64-apple-darwin20
> Running under: macOS 15.0.1
>
> Matrix products: default
> BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> time zone: America/Phoenix
> tzcode source: internal
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] vioplot_0.5.0 zoo_1.8-12 sm_2.2-6.0 scales_1.3.0 rphylopic_1.4.0.9000 raster_3.6-26
> [7] sp_2.1-4 nlme_3.1-165 MuMIn_1.48.4 lattice_0.22-6 phylolm_2.6.2 geiger_2.0.11
> [13] phytools_2.3-0 maps_3.4.2 ape_5.8 AICcmodavg_2.3-3 rmarkdown_2.27
>
> loaded via a namespace (and not attached):
> [1] tidyselect_1.2.1 subplex_1.9 farver_2.1.2 dplyr_1.1.4 optimParallel_1.0-2
> [6] fastmap_1.2.0 combinat_0.0-8 XML_3.99-0.17 digest_0.6.36 lifecycle_1.0.4
> [11] rsvg_2.6.0 survival_3.6-4 terra_1.7-78 magrittr_2.0.3 compiler_4.4.1
> [16] rlang_1.1.4 sass_0.4.9 tools_4.4.1 igraph_2.0.3 utf8_1.2.4
> [21] yaml_2.3.8 knitr_1.47 phangorn_2.11.1 clusterGeneration_1.3.8 grImport2_0.3-1
> [26] mnormt_2.1.1 scatterplot3d_0.3-44 curl_5.2.1 expm_0.999-9 numDeriv_2016.8-1.1
> [31] grid_4.4.1 stats4_4.4.1 fansi_1.0.6 unmarked_1.4.1 xtable_1.8-4
> [36] colorspace_2.1-0 future_1.34.0 ggplot2_3.5.1 globals_0.16.3 iterators_1.0.14
> [41] MASS_7.3-60.2 cli_3.6.3 mvtnorm_1.2-5 generics_0.1.3 future.apply_1.11.2
> [46] httr_1.4.7 pbapply_1.7-2 cachem_1.1.0 splines_4.4.1 parallel_4.4.1
> [51] base64enc_0.1-3 vctrs_0.6.5 Matrix_1.7-0 jsonlite_1.8.8 VGAM_1.1-11
> [56] listenv_0.9.1 jpeg_0.1-10 foreach_1.5.2 jquerylib_0.1.4 glue_1.7.0
> [61] parallelly_1.38.0 codetools_0.2-20 DEoptim_2.2-8 gtable_0.3.5 quadprog_1.5-8
> [66] munsell_0.5.1 tibble_3.2.1 pillar_1.9.0 htmltools_0.5.8.1 deSolve_1.40
> [71] R6_2.5.1 doParallel_1.0.17 evaluate_0.24.0 highr_0.11 png_0.1-8
> [76] bslib_0.7.0 Rcpp_1.0.12 fastmatch_1.1-4 coda_0.19-4.1 xfun_0.45
> [81] pkgconfig_2.0.3Beginning of the analyses
Figure 1.
#png("tree.png", height = 7, width = 7, units = "in", res = 360)
col.br <- setNames(c("purple", "orange"), c("Archaea", "Bacteria"))
plotTree(tree.vol, ftype = "i", lwd = 3, mar = c(3.5, 1, 1, 3))
par(new = TRUE, col = "transparent")
painted <- paintSubTree(tree.vol, 52, "Archaea" ,"0")
plotSimmap(painted, col.br, ftype = "i", lwd = 3, mar = c(3.5, 1, 1, 3))
par(new = TRUE, col = "transparent")
painted <- paintSubTree(tree.vol, 31, "Bacteria")
plotSimmap(painted, col.br, ftype = "i", lwd = 3, mar = c(3.5, 1, 1, 3))
par(new = TRUE, col = "black")
legend("bottomleft", legend = c("Archaea", "Bacteria"), lwd = 3, col = col.br, bty = "n")
axisPhylo(1, line = -0.1)
mtext("Time (mya)", side = 1, line = 2, at = 2000)
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x2 <- runif(100, obj$x.lim[2] + 10, obj$x.lim[2] + 50)
spp <- gsub("(_).*","", tree.vol$tip.label)[-c(3, 6, 9, 10, 12, 14, 16, 18, 21, 24, 28)]
spp[7] <- "Cyanobacteria"
spp[8] <- "Mycoplasma_genitalium"
spp[10] <- "Pleurocapsa fuliginosa"
col.fill <- c(rep("orange", 13), rep("purple", 7))
col.bor <- c(rep("black", 13), rep("red", 7))
idx <- 1
for(i in spp){
x2
uuid <- get_uuid(name = i, n = 1)
img <- get_phylopic(uuid = uuid)
nodes <- sapply(i, grep, x = tree.vol$tip.label)
for(j in nodes){
add_phylopic_base(img = img, x = sample(x2, 1), y = j, ysize = 1, color = col.bor[idx], fill = col.fill[idx])
}
idx = idx + 1
}
uuid <- get_uuid(name = "Chroococcus turgidus", n = 1)
img <- get_phylopic(uuid = uuid)
add_phylopic_base(img = img, x = 7400, y = 10, ysize = 1, color = "black", fill = "orange")
uuid <- get_uuid(name = "Pleurocapsa fuliginosa", n = 1)
img <- get_phylopic(uuid = uuid)
add_phylopic_base(img = img, x = 7400, y = 17, ysize = 1, color = "black", fill = "orange")
uuid <- get_uuid(name = "Fimbriimonas ginsengisoli", n = 1)
img <- get_phylopic(uuid = uuid)
add_phylopic_base(img = img, x = 7400, y = 11.5, ysize = 1, color = alpha("black", 0.1), fill = "orange")
#dev.off()Figure 4
tmp.lower.dat <- tmp.lower.dat[tmp.lower.dat$d1_up < 4, ] ## Removing potential outliers
tmp.op.dat <- tmp.op.dat[tmp.op.dat$d1_up < 4, ] ## Removing potential outliers
tmp.upper.dat <- tmp.upper.dat[tmp.upper.dat$d1_up < 4, ] ## Removing potential outliers
#png("figure4.png", height = 7, width = 7, units = "in", res = 360)
layout(matrix(c(1, 1, 2, 2, 3, 3,
1, 1, 2, 2, 3, 3,
4, 4, 5, 5, 6, 6,
4, 4, 5, 5, 6, 6), nrow = 4, ncol = 6, byrow = TRUE))
## lower
tmp.lower.dat <- tmp.lower.dat[tmp.lower.dat$d1_up < 4, ] ## Removing potential outliers
model6.1 <- lm(d1_up ~ log10(tmp.op), data = tmp.lower.dat)
##summary(model6.1)
## IC
SSX <- sum(round((log10(tmp.lower.dat$tmp.op) - mean(log10(tmp.lower.dat$tmp.op)))^2), 2)
s2 <- var(tmp.lower.dat$d1_up)
n <- length(tmp.lower.dat$d1_up)
x <- seq(min(log10(tmp.lower.dat$tmp.op)), max(log10(tmp.lower.dat$tmp.op)), length = 30)
m.x <- mean(round(log10(tmp.lower.dat$tmp.op), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model6.1)[1] + coef(model6.1)[2]*x) + ic.s
lower.i <- (coef(model6.1)[1] + coef(model6.1)[2]*x) + ic.i
par(mar = c(6.4, 4, 2, 0), mgp = c(2.8, 1, 0))
plot(d1_up ~ log10(tmp.op), data = tmp.lower.dat, ylab = expression(paste("Cell diameter ", log[10], sep = " ")(mu*m)), xlab = expression(paste("Lower temperature")~log[10]~("\u00B0C")), las = 1, pch = 21, cex = 1.2, type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
cols <- setNames(c("purple", "orange"), levels(as.factor(tmp.lower.dat$kingdom)))
plot(d1_up ~ log10(tmp.op), data = tmp.lower.dat, xlab = "", ylab = "", las = 1, pch = 21, col = cols[tmp.lower.dat$kingdom], bg = cols[tmp.lower.dat$kingdom], cex = 0.8, axes = FALSE)
#lines(x = x, y = (coef(model6.1)[1] + coef(model6.1)[2]*x), lwd = 2, col = "black")
#polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
legend("topleft", legend = levels(as.factor((tmp.lower.dat$kingdom))), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 0.8)
legend("topright", legend = paste("n = ", length(tmp.lower.dat$species)), bty = "n", cex = 0.8)
mtext("A", side = 2, line = 2.6, at = 3.8, las = 1, font = 2)
## Cell size and temp opt
tmp.op.dat <- tmp.op.dat[tmp.op.dat$d1_up < 4, ] ## Removing potential outliers
#model6 <- gls(d1_up ~ log10(tmp.op), correlation = corBrownian(phy = tree.tmp, form = ~species), data = tmp.dat, method = "ML")
#model6 <- lm(d1_up ~ log10(tmp.op), data = tmp.dat)
model6 <- lm(d1_up ~ log10(tmp.op), data = tmp.op.dat)
#summary(model6)
## IC
SSX <- sum(round((log10(tmp.op.dat$tmp.op) - mean(log10(tmp.op.dat$tmp.op)))^2), 2)
s2 <- var(tmp.op.dat$d1_up)
n <- length(tmp.op.dat$d1_up)
x <- seq(min(log10(tmp.op.dat$tmp.op)), max(log10(tmp.op.dat$tmp.op)), length = 30)
m.x <- mean(round(log10(tmp.op.dat$tmp.op), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model6)[1] + coef(model6)[2]*x) + ic.s
lower.i <- (coef(model6)[1] + coef(model6)[2]*x) + ic.i
par(mar = c(6.4, 2.3, 2, 0.1))
plot(d1_up ~ log10(tmp.op), data = tmp.op.dat, ylab = " ", xlab = expression(paste("Optimum temperature")~log[10]~("\u00B0C")), las = 1, pch = 21, cex = 1.2, type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
cols <- setNames(c("purple", "orange"), levels(as.factor(tmp.op.dat$kingdom)))
plot(d1_up ~ log10(tmp.op), data = tmp.op.dat, xlab = "", ylab = "", las = 1, pch = 21, col = cols[tmp.op.dat$kingdom], bg = cols[tmp.op.dat$kingdom], cex = 0.8, axes = FALSE)
#lines(x = x, y = (coef(model6)[1] + coef(model6)[2]*x), lwd = 2, col = "black")
#polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
legend("topleft", legend = levels(as.factor((tmp.op.dat$kingdom))), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 0.8)
legend(x = 0.7, y = 3.07, legend = paste("n = ", length(tmp.op.dat$species)), bty = "n", cex = 0.8)
## upper
tmp.upper.dat <- tmp.upper.dat[tmp.upper.dat$d1_up < 4, ] ## Removing potential outliers
model6.3 <- lm(d1_up ~ log10(tmp.op), data = tmp.upper.dat)
#summary(model6.3)
## Model selection procedure based on AIC criterion
AICc(model6.1, model6, model6.3)> df AICc
> model6.1 3 751.0115
> model6 3 405.3537
> model6.3 3 519.2041## IC
SSX <- sum(round((log10(tmp.upper.dat$tmp.op) - mean(log10(tmp.upper.dat$tmp.op)))^2), 2)
s2 <- var(tmp.upper.dat$d1_up)
n <- length(tmp.upper.dat$d1_up)
x <- seq(min(log10(tmp.upper.dat$tmp.op)), max(log10(tmp.upper.dat$tmp.op)), length = 30)
m.x <- mean(round(log10(tmp.upper.dat$tmp.op), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model6.3)[1] + coef(model6.3)[2]*x) + ic.s
lower.i <- (coef(model6.3)[1] + coef(model6.3)[2]*x) + ic.i
par(mar = c(6.4, 2.3, 2, 0.2))
plot(d1_up ~ log10(tmp.op), data = tmp.upper.dat, ylab = " ", xlab = expression(paste("Upper temperature")~log[10]~("\u00B0C")), las = 1, pch = 21, cex = 1.2, type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
cols <- setNames(c("purple", "orange"), levels(as.factor(tmp.upper.dat$kingdom)))
plot(d1_up ~ log10(tmp.op), data = tmp.upper.dat, xlab = "", ylab = "", las = 1, pch = 21, col = cols[tmp.upper.dat$kingdom], bg = cols[tmp.upper.dat$kingdom], cex = 0.8, axes = FALSE)
lines(x = x, y = (coef(model6.3)[1] + coef(model6.3)[2]*x), lwd = 2, col = "black")
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
legend("topleft", legend = levels(as.factor((tmp.upper.dat$kingdom))), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 0.8)
legend(x = 0.6, y = 3, legend = paste("n = ", length(tmp.upper.dat$species)), bty = "n", cex = 0.8)
legend(x = 0.6, y = 2.85, legend = expression(p==0.002), bty = "n", cex = 0.8)
legend(x = 0.6, y = 2.7, legend = expression(R^2==0.02), bty = "n", cex = 0.8)
## GROWTH
## lower
tmp.lower.dat.growth$r <- log(2)/tmp.lower.dat.growth$doubling_h
model7.1 <- lm(log10(r) ~ log10(tmp.op), data = tmp.lower.dat.growth)
## IC
SSX <- sum(round((log10(tmp.lower.dat.growth$tmp.op) - mean(log10(tmp.lower.dat.growth$tmp.op)))^2), 2)
s2 <- var(log10(tmp.lower.dat.growth$r))
n <- length(log10(tmp.lower.dat.growth$r))
x <- seq(min(log10(tmp.lower.dat.growth$tmp.op)), max(log10(tmp.lower.dat.growth$tmp.op)), length = 30)
m.x <- mean(round(log10(tmp.lower.dat.growth$tmp.op), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model7.1)[1] + coef(model7.1)[2]*x) + ic.s
lower.i <- (coef(model7.1)[1] + coef(model7.1)[2]*x) + ic.i
par(mar = c(6.4, 4, 2, 0), mgp = c(2.7, 1, 0))
plot(log10(r) ~ log10(tmp.op), data = tmp.lower.dat.growth, ylab = expression(paste('Intrinsic growth rate')~italic(r)~log[10]~(h^-1)), xlab = expression(paste("Lower temperature")~log[10]~("\u00B0C")), las = 1, pch = 21, bg = alpha("black", 0.3), cex = 1.2, type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
cols3 <- c("purple", "orange")[as.numeric(as.factor(tmp.lower.dat.growth$kingdom))]
plot(log10(r) ~ log10(tmp.op), data = tmp.lower.dat.growth, xlab = "", ylab = "", las = 1, pch = 21, col = cols3, bg = cols3, cex = 0.8, axes = FALSE)
lines(x = x, y = (coef(model7.1)[1] + coef(model7.1)[2]*x), lty = 2, lwd = 2, col = "black")
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
legend("bottomright", legend = levels(as.factor((tmp.lower.dat.growth$kingdom))), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 0.8)
mtext("B", side = 2, line = 2.6, at = 0.7, las = 1, font = 2)
legend(x = 1.13, y = 0.35, legend = paste("n = ", length(tmp.lower.dat.growth$species)), bty = "n", cex = 0.8)
legend(x = 1.13, y = 0.2, legend = expression(p<0.001), bty = "n", cex = 0.8)
legend(x = 1.13, y = 0.07, legend = expression(R^2==0.42), bty = "n", cex = 0.8)
## Doubling and optimum temp
#model7 <- gls(log(doubling_h) ~ log(tmp.op), correlation = corBrownian(phy = tree.tmp2, form = ~species), data = tmp.dat2, method = "ML")
tmp.op.dat.growth$r <- log(2)/tmp.op.dat.growth$doubling_h
model7 <- lm(log10(r) ~ log10(tmp.op), data = tmp.op.dat.growth)
#summary(model7)
## IC
SSX <- sum(round((log10(tmp.op.dat.growth$tmp.op) - mean(log10(tmp.op.dat.growth$tmp.op)))^2), 2)
s2 <- var(log10(tmp.op.dat.growth$r))
n <- length(log10(tmp.op.dat.growth$r))
x <- seq(min(log10(tmp.op.dat.growth$tmp.op)), max(log10(tmp.op.dat.growth$tmp.op)), length = 30)
m.x <- mean(round(log10(tmp.op.dat.growth$tmp.op), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model7)[1] + coef(model7)[2]*x) + ic.s
lower.i <- (coef(model7)[1] + coef(model7)[2]*x) + ic.i
par(mar = c(6.4, 2.3, 2, 0.1))
plot(log10(r) ~ log10(tmp.op), data = tmp.op.dat.growth, ylab = expression(paste("Intrinsic growth rate")~r~log[10]~(h^-1)), xlab = expression(paste("Optimum temperature")~log[10]~("\u00B0C")), las = 1, pch = 21, bg = alpha("black", 0.3), cex = 1.2, type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
cols3 <- c("purple", "orange")[as.numeric(as.factor(tmp.op.dat.growth$kingdom))]
plot(log10(r) ~ log10(tmp.op), data = tmp.op.dat.growth, xlab = "", ylab = "", las = 1, pch = 21, col = cols3, bg = cols3, cex = 0.8, axes = FALSE)
lines(x = x, y = (coef(model7)[1] + coef(model7)[2]*x), lty = 2, lwd = 2, col = "black")
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
legend("topleft", legend = levels(as.factor((tmp.op.dat.growth$kingdom))), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 0.8)
legend(x = 0.92, y = -0.01, legend = paste("n = ", length(tmp.op.dat.growth$species)), bty = "n", cex = 0.8)
legend(x = 0.92, y = -0.12, legend = expression(p==0.02), bty = "n", cex = 0.8)
legend(x = 0.92, y = -0.21, legend = expression(R^2==0.21), bty = "n", cex = 0.8)
## upper
tmp.upper.dat.growth$r <- log(2)/tmp.upper.dat.growth$doubling_h
model7.2 <- lm(log10(r) ~ log10(tmp.op), data = tmp.upper.dat.growth)
## IC
SSX <- sum(round((log10(tmp.upper.dat.growth$tmp.op) - mean(log10(tmp.upper.dat.growth$tmp.op)))^2), 2)
s2 <- var(log10(tmp.upper.dat.growth$r))
n <- length(log10(tmp.upper.dat.growth$r))
x <- seq(min(log10(tmp.upper.dat.growth$tmp.op)), max(log10(tmp.upper.dat.growth$tmp.op)), length = 30)
m.x <- mean(round(log10(tmp.upper.dat.growth$tmp.op), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(model7.2)[1] + coef(model7.2)[2]*x) + ic.s
lower.i <- (coef(model7.2)[1] + coef(model7.2)[2]*x) + ic.i
par(mar = c(6.4, 2.3, 2, 0.2))
plot(log10(r) ~ log10(tmp.op), data = tmp.upper.dat.growth, ylab = expression(paste("Intrinsic growth rate")~r~log[10]~(h^-1)), xlab = expression(paste("Upper temperature")~log[10]~("\u00B0C")), las = 1, pch = 21, bg = alpha("black", 0.3), cex = 1.2, type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
cols3 <- c("purple", "orange")[as.numeric(as.factor(tmp.upper.dat.growth$kingdom))]
plot(log10(r) ~ log10(tmp.op), data = tmp.upper.dat.growth, xlab = "", ylab = "", las = 1, pch = 21, col = cols3, bg = cols3, cex = 0.8, axes = FALSE)
lines(x = x, y = (coef(model7.2)[1] + coef(model7.2)[2]*x), lty = 2, lwd = 2, col = "black")
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
legend("topleft", legend = levels(as.factor((tmp.upper.dat.growth$kingdom))), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 0.8)
legend(x = 0.6, y = 0.16, legend = paste("n = ", length(tmp.upper.dat.growth$species)), bty = "n", cex = 0.8)
legend(x = 0.6, y = 0.05, legend = expression(p<0.001), bty = "n", cex = 0.8)
legend(x = 0.6, y = -0.04, legend = expression(R^2==0.22), bty = "n", cex = 0.8)
#dev.off()
## int.mod <- lm(log10(r) ~ log10(tmp.op)*d1_up, data = tmp.op.dat.growth) ## This is the model that describes the interaction in question
## summary(int.mod)
##AICc(model6, model6.1, model6.3, model7, model7.1, int.mod)
##anova(model7, int.mod)Model diagnosis
model6.3 <- lm(d1_up ~ log10(tmp.op), data = tmp.upper.dat)
summary(model6.3)>
> Call:
> lm(formula = d1_up ~ log10(tmp.op), data = tmp.upper.dat)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.80822 -0.33283 -0.13363 0.09665 2.31893
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.04328 0.29301 0.148 0.88267
> log10(tmp.op) 0.61525 0.19930 3.087 0.00219 **
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.5281 on 328 degrees of freedom
> Multiple R-squared: 0.02823, Adjusted R-squared: 0.02527
> F-statistic: 9.53 on 1 and 328 DF, p-value: 0.002194layout(matrix(c(1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
par(las = 1)
plot(model6.3)
Figure 2
v.dat$r <- log(2)/v.dat$doubling_h
mod <- lm(log10(r) ~ log10(volume), data = v.dat)
#summary(mod)
pg.mod <- gls(log10(r) ~ log10(volume), correlation = corBrownian(phy = tree.vol, form = ~species), data = vol.dat, method = "ML")
#summary(pg.mod)
AICc(mod, pg.mod)> df AICc
> mod 3 53.55435
> pg.mod 3 18.57799#png("figure2.png", height = 7, width = 7, units = "in", res = 360)
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## IC
SSX <- sum(round((log10(vol.dat$volume) - mean(log10(vol.dat$volume)))^2), 2)
s2 <- var(log10(vol.dat$r))
n <- length(vol.dat$r)
x <- seq(min(log10(vol.dat$volume)), max(log10(vol.dat$volume)), length = length(vol.dat$species))
m.x <- mean(round(log(vol.dat$volume), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(pg.mod)[1] + coef(pg.mod)[2]*x) + ic.s
lower.i <- (coef(pg.mod)[1] + coef(pg.mod)[2]*x) + ic.i
par(mgp = c(2.5, 1, 0))
cols <- setNames(c("purple", "orange"), levels(as.factor(vol.dat$kingdom)))
vol.dat$kingdom <- as.factor(vol.dat$kingdom)
plot(log10(r) ~ log10(volume), data = vol.dat, type = "n", pch = 21, las = 1, ylab = expression(paste('Intrinsic growth rate')~italic(r)~log[10]~(h^-1)), xlab = expression(paste("Cell volume ", log[10], sep = " ")(mu*m^3)))
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
plot(log10(r) ~ log10(volume), data = vol.dat, type = "p", pch = 21, col = cols[vol.dat$kingdom], bg = cols[vol.dat$kingdom], las = 1, ylab = "", xlab = "", axes = FALSE)
lines(x, y = (coef(pg.mod)[1] + coef(pg.mod)[2]*x), lwd = 2)
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
legend("topleft", legend = unique(vol.dat$kingdom), pch = 16, col = cols, bg = cols, bty = "n")
legend(x = -2.7, y = 0, legend = paste('n = ', length(vol.dat$r), sep = ''), bty = "n")
legend(x = -2.7, y = -0.13, legend = expression(p<0.001), bty = 'n')
mtext("A", side = 2, at = 0.7, line = 1, las = 1, font = 2)
## Doubling and optimum temp
mod.fg1 <- lm(log10(r) ~ d1_up, data = tmp.op.dat.growth)
#summary(mod.fg1)
check <- name.check(tree, tmp.op.dat.growth)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
pg2.tree <- drop.tip(tree, rm_phy)
pg.mod2 <- gls(log10(r) ~ d1_up, correlation = corBrownian(phy = pg2.tree , form = ~species), data = tmp.op.dat.growth, method = "ML")
summary(pg.mod2)> Generalized least squares fit by maximum likelihood
> Model: log10(r) ~ d1_up
> Data: tmp.op.dat.growth
> AIC BIC logLik
> 138.4652 143.5318 -66.2326
>
> Correlation Structure: corBrownian
> Formula: ~species
> Parameter estimate(s):
> numeric(0)
>
> Coefficients:
> Value Std.Error t-value p-value
> (Intercept) -1.5787613 0.8887161 -1.776452 0.0837
> d1_up 0.7831775 0.2289400 3.420885 0.0015
>
> Correlation:
> (Intr)
> d1_up -0.255
>
> Standardized residuals:
> Min Q1 Med Q3 Max
> -1.11225248 -0.27014620 -0.08146687 0.04575048 0.29891957
>
> Residual standard error: 2.361099
> Degrees of freedom: 40 total; 38 residualAICc(mod.fg1, pg.mod2)> df AICc
> mod.fg1 3 78.15646
> pg.mod2 3 139.13187#xtable(summary(mod.fg1)$coefficients, digits = 3)
par(mgp = c(2.5, 1, 0))
plot(log10(r) ~ d1_up, data = tmp.op.dat.growth, ylab = expression(paste("Intrinsic growth rate")~italic(r)~log[10]~(h^-1)), xlab = expression(paste("Cell diameter")~log[10]~(mu*m)), las = 1, pch = 21, bg = alpha("black", 0.3), cex = 1.2, type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
cols3 <- c("purple", "orange")[as.numeric(as.factor(tmp.op.dat.growth$kingdom))]
plot(log10(r) ~ d1_up, data = tmp.op.dat.growth, xlab = "", ylab = "", las = 1, pch = 21, col = cols3, bg = cols3, cex = 1, axes = FALSE)
legend("topleft", legend = levels(as.factor((tmp.op.dat.growth$kingdom))), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 1)
mtext("B", side = 2, at = 0.48, line = 1, las = 1, font = 2)
#dev.off()Model diagnosis
mod.fg1 <- lm(log10(doubling_h) ~ d1_up, data = tmp.op.dat.growth)
summary(mod.fg1)>
> Call:
> lm(formula = log10(doubling_h) ~ d1_up, data = tmp.op.dat.growth)
>
> Residuals:
> Min 1Q Median 3Q Max
> -0.93943 -0.42732 0.02871 0.46205 1.36485
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1.1456 0.1689 6.782 4.86e-08 ***
> d1_up -0.1691 0.1392 -1.214 0.232
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.6067 on 38 degrees of freedom
> Multiple R-squared: 0.03737, Adjusted R-squared: 0.01203
> F-statistic: 1.475 on 1 and 38 DF, p-value: 0.2321layout(matrix(c(1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
par(las = 1)
plot(mod.fg1)
Figure 3
## Translation machinery
tRNA <- aggregate(spp.tRNA$tRNA_genes, by = list(spp.tRNA$species), mean, na.action = na.rm)
rRNA <- aggregate(spp.rRNA$rRNA16S_genes, by = list(spp.rRNA$species), mean, na.action = na.rm)
cell.vol <- aggregate(vol$volume, by = list(vol$species), mean)
d1_up <- aggregate(spp.d1_up$d1_up, by = list(spp.d1_up$species), mean, na.action = na.rm)
doubling <- aggregate(spp.doubling$doubling_h, by = list(spp.doubling$species), mean, na.action = na.rm)
#dim(tRNA)
names(tRNA) <- c("species", "tRNA")
#dim(rRNA)
names(rRNA) <- c("species", "rRNA")
#dim(cell.vol)
names(cell.vol) <- c("species", "volume")
#dim(d1_up)
names(d1_up) <- c("species", "d1_up")
#dim(doubling)
names(doubling) <- c("species", "doubling_h")
genes <- merge(rRNA, tRNA, by = "species")
tran <- merge(genes, cell.vol, by = "species")
tran2 <- merge(genes, d1_up, by = "species")
tran3 <- merge(genes, doubling, by = "species")
obj <- rep()
for(i in tran$species){
kingdom <- unique(data$superkingdom[data$species == i])
obj <- c(obj, kingdom)
}
tran$kingdom <- obj
#head(tran)
tran$species <- gsub("[[:punct:]]", "", tran$species)
tran$species <- gsub(" ", "_", tran$species)
#head(tran)
tran <- tran[!tran$species == "Sphingopyxis_alaskensis", ] ## possible outlier
rownames(tran) <- tran$species
tran.tree <- read.tree("tran.spp.nwk")
#tran.tree <- force.ultrametric(tran.tree)
check <- name.check(tran.tree, tran)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
tree.tran <- drop.tip(tran.tree, rm_phy)
tran.dat <- subset(tran, subset = tran$species %in% tree.tran$tip, select = names(tran))
#name.check(tree.tran, tran.dat)
#str(tran.dat)
obj <- rep()
for(i in tran2$species){
kingdom <- unique(data$superkingdom[data$species == i])
obj <- c(obj, kingdom)
}
tran2$kingdom <- obj
#head(tran)
tran2$species <- gsub("[[:punct:]]", "", tran2$species)
tran2$species <- gsub(" ", "_", tran2$species)
rownames(tran2) <- tran2$species
#head(tran2)
check <- name.check(tree, tran2)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
tree.tran2 <- drop.tip(tree, rm_phy)
tran.dat2 <- subset(tran2, subset = tran2$species %in% tree$tip, select = names(tran2))
#name.check(tree.tran2, tran.dat2)
#str(tran.dat2)
obj <- rep()
for(i in tran3$species){
kingdom <- unique(data$superkingdom[data$species == i])
obj <- c(obj, kingdom)
}
tran3$kingdom <- obj
#head(tran)
tran3$species <- gsub("[[:punct:]]", "", tran3$species)
tran3$species <- gsub(" ", "_", tran3$species)
rownames(tran3) <- tran3$species
tran3$r <- log(2)/tran3$doubling_h
#head(tran3)
check <- name.check(tree, tran3)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
tree.tran3 <- drop.tip(tree, rm_phy)
tran.dat3 <- subset(tran3, subset = tran3$species %in% tree.tran3$tip, select = names(tran3))
##str(tran.dat3)
mod.gr1 <- lm(log10(rRNA) ~ log10(r), data = tran3)
#summary(mod.gr1)
mod.gr2 <- lm(log10(tRNA) ~ log10(r), data = tran3)
#summary(mod.gr2)
int.mod <- lm(log10(tRNA) ~ log10(r)*log10(rRNA), data = tran3)
#summary(int.mod)
int.mod2 <- lm(log10(rRNA) ~ log10(r)*log10(tRNA), data = tran3)
#summary(int.mod2)
anova(mod.gr2, int.mod)> Analysis of Variance Table
>
> Model 1: log10(tRNA) ~ log10(r)
> Model 2: log10(tRNA) ~ log10(r) * log10(rRNA)
> Res.Df RSS Df Sum of Sq F Pr(>F)
> 1 411 5.7439
> 2 409 2.0998 2 3.6441 354.91 < 2.2e-16 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICc(mod.gr2, int.mod, int.mod2)> df AICc
> mod.gr2 3 -587.6016
> int.mod 5 -999.1221
> int.mod2 5 -187.7779#xtable(summary(int.mod)$coefficients, digits = 3)
## GROWTH
#png("figure3.png", height = 7, width = 7, units = "in", res = 360)
SSX <- sum(round((log10(tran3$r) - mean(log10(tran3$r)))^2), 2)
s2 <- var(log10(tran3$rRNA))
n <- length(tran3$rRNA)
x <- seq(min(log10(tran3$r)), max(log10(tran3$r)), length = length(tran3$species))
m.x <- mean(round(log(tran3$r), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(mod.gr1)[1] + coef(mod.gr1)[2]*x) + ic.s
lower.i <- (coef(mod.gr1)[1] + coef(mod.gr1)[2]*x) + ic.i
cols2 <- setNames(c("purple", "orange"), levels(as.factor(tran3$kingdom)))
par(mar = c(5, 4, 1.5, 3.3), mgp = c(2.7, 1, 0))
plot(log10(rRNA) ~ log10(r), data = tran3, type = "n", pch = 16, las = 1, ylab = expression(log[10]~rRNA~genes), xlab = expression(paste("Intrinsic growth rate")~italic(r)~log[10]~(h^-1)))
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
plot(log10(rRNA) ~ log10(r), data = tran3, type = "p", pch = 16, col = cols2[tran3$kingdom], bg = cols2[tran3$kingdom], las = 1, axes = FALSE, xaxt = "n", ylab = "", xlab = "")
lines(x = x, y = (coef(mod.gr1)[1] + coef(mod.gr1)[2]*x), lwd = 2)
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
SSX <- sum(round((log10(tran3$r) - mean(log10(tran3$r)))^2), 2)
s2 <- var(log10(tran3$tRNA))
n <- length(tran3$tRNA)
x <- seq(min(log10(tran3$r)), max(log10(tran3$r)), length = length(tran3$species))
m.x <- mean(round(log(tran3$rRNA), 1))
se <- sqrt(s2*((1/n) + (((x - m.x)^2)/SSX)))
is <- qt(0.975, df = n - 2)
ii <- qt(0.025, df = n - 2)
ic.s <- se*is
ic.i <- se*ii
upper.i <- (coef(mod.gr2)[1] + coef(mod.gr2)[2]*x) + ic.s
lower.i <- (coef(mod.gr2)[1] + coef(mod.gr2)[2]*x) + ic.i
cols2 <- setNames(c("purple", "orange"), levels(as.factor(tran3$kingdom)))
#plot(log10(tRNA) ~ log10(doubling_h), data = tran3, type = "n", pch = 16, las = 1, ylab = expression(log[10]~tRNA~genes), xlab = expression(paste("Doubling ", log[10], sep = " ")(h)))
#grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
plot(log10(tRNA) ~ log10(r), data = tran3, type = "p", pch = 1, col = cols2[tran3$kingdom], bg = cols2[tran3$kingdom], las = 1, axes = FALSE, xaxt = "n", ylab = "", xlab = "")
lines(x = x, y = (coef(mod.gr2)[1] + coef(mod.gr2)[2]*x), lwd = 2, lty = 2)
polygon(c(rev(x), x), c(rev(lower.i), upper.i), border = FALSE, col = alpha("gold", 0.3))
axis(side = 4, at = pretty(range(log10(tran3$tRNA))), las = 1)
mtext(expression(log[10]~tRNA~genes), side = 4, line = 2.3)
legend("topleft", legend = c("Archaea", "Bacteria"), pch = 16, col = c("purple", "orange"), bg = c("purple", "orange"), bty = "n", cex = 1)
legend(x = -2.9, y = 2.19, legend = "tRNA genes", lty = 2, lwd = 2, bty = "n")
legend(x = -2.94, y = 2.259, legend = c(" ", " "), pch = 1, col = c("purple", "orange"), bty = "n", cex = 1)
#legend(x = -3, y = 2.15, legend = expression(p==0.04), bty = "n")
#legend(x = -3, y = 2.12, legend = expression(R^2==0.17), bty = "n")
#dev.off()Model diagnosis
layout(matrix(c(1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
par(las = 1)
plot(int.mod)
Figure 5
obj <- rep()
for(i in d1_up$species){
kingdom <- unique(data$superkingdom[data$species == i])
obj <- c(obj, kingdom)
}
d1_up$kingdom <- obj
head(d1_up)> species d1_up kingdom
> 1 [Clostridium] aldenense 1.1 Bacteria
> 2 [Clostridium] caenicola 0.6 Bacteria
> 3 [Clostridium] fimetarium 0.6 Bacteria
> 4 [Clostridium] lavalense 1.5 Bacteria
> 5 [Clostridium] paradoxum 1.1 Bacteria
> 6 [Clostridium] polysaccharolyticum 1.1 Bacteriad1_up$species <- gsub("[[:punct:]]", "", d1_up$species)
d1_up$species <- gsub(" ", "_", d1_up$species)
rownames(d1_up) <- d1_up$species
dim(d1_up)> [1] 1603 3head(d1_up)> species d1_up kingdom
> Clostridium_aldenense Clostridium_aldenense 1.1 Bacteria
> Clostridium_caenicola Clostridium_caenicola 0.6 Bacteria
> Clostridium_fimetarium Clostridium_fimetarium 0.6 Bacteria
> Clostridium_lavalense Clostridium_lavalense 1.5 Bacteria
> Clostridium_paradoxum Clostridium_paradoxum 1.1 Bacteria
> Clostridium_polysaccharolyticum Clostridium_polysaccharolyticum 1.1 Bacteriad1_up <- d1_up[d1_up$d1_up < 6, ]
size <- lm(d1_up ~ kingdom, data = d1_up)
summary(size)>
> Call:
> lm(formula = d1_up ~ kingdom, data = d1_up)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1.1004 -0.2783 -0.0783 0.1217 4.6217
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1.25044 0.04059 30.803 <2e-16 ***
> kingdomBacteria -0.37217 0.04339 -8.577 <2e-16 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.5727 on 1594 degrees of freedom
> Multiple R-squared: 0.04412, Adjusted R-squared: 0.04352
> F-statistic: 73.57 on 1 and 1594 DF, p-value: < 2.2e-16obj <- rep()
for(i in doubling$species){
kingdom <- unique(data$superkingdom[data$species == i])
obj <- c(obj, kingdom)
}
doubling$kingdom <- obj
head(doubling)> species doubling_h kingdom
> 1 [Bacillus] caldotenax 0.24 Bacteria
> 2 [Butyribacterium] methylotrophicum 20.00 Bacteria
> 3 [Clostridium] alkalicellulosi 14.00 Bacteria
> 4 [Clostridium] paradoxum 0.67 Bacteria
> 5 [Clostridium] stercorarium 8.60 Bacteria
> 6 [Clostridium] termitidis 9.62 Bacteriadoubling$species <- gsub("[[:punct:]]", "", doubling$species)
doubling$species <- gsub(" ", "_", doubling$species)
rownames(doubling) <- doubling$species
doubling$r <- log(2)/doubling$doubling_h
dim(doubling)> [1] 928 4head(doubling)> species doubling_h kingdom r
> Bacillus_caldotenax Bacillus_caldotenax 0.24 Bacteria 2.88811325
> Butyribacterium_methylotrophicum Butyribacterium_methylotrophicum 20.00 Bacteria 0.03465736
> Clostridium_alkalicellulosi Clostridium_alkalicellulosi 14.00 Bacteria 0.04951051
> Clostridium_paradoxum Clostridium_paradoxum 0.67 Bacteria 1.03454803
> Clostridium_stercorarium Clostridium_stercorarium 8.60 Bacteria 0.08059851
> Clostridium_termitidis Clostridium_termitidis 9.62 Bacteria 0.07205272growth <- lm(log10(r) ~ kingdom, data = doubling)
summary(growth)>
> Call:
> lm(formula = log10(r) ~ kingdom, data = doubling)
>
> Residuals:
> Min 1Q Median 3Q Max
> -2.10675 -0.46236 0.05262 0.45073 1.37954
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) -0.90044 0.04201 -21.435 < 2e-16 ***
> kingdomBacteria 0.21152 0.04826 4.383 1.31e-05 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.6301 on 926 degrees of freedom
> Multiple R-squared: 0.02032, Adjusted R-squared: 0.01926
> F-statistic: 19.21 on 1 and 926 DF, p-value: 1.307e-05##png("figure5.png", height = 7, width = 7, units = "in", res = 360)
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
0, 0, 0, 0), nrow = 4, ncol = 4, byrow = TRUE))
## Basic boxplot
par(mgp = c(2.5, 1, 0))
vioplot(log10(r) ~ kingdom, data = doubling, border = NA, method = "jitter", side = "right", ylab = expression(paste("Intrinsic growth rate")~italic(r)~log[10]~(h^-1)), xlab = "Kingdom", col = "white", las = 1)
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
vioplot(log10(r) ~ kingdom, data = doubling, border = NA, method = "jitter", side = "right", ylab = "", xlab = "", col = c(alpha("purple", 0.2), alpha("orange", 0.2)), las = 1)
segments(x0 = 1, y0 = mean(log10(doubling$r)[doubling$kingdom == "Archaea"]), x1 = 1.4, y1 = mean(log10(doubling$r)[doubling$kingdom == "Archaea"]), lwd = 2, lty = 2, col = "black")
text(x = 1.45, y = mean(log10(doubling$r)[doubling$kingdom == "Archaea"]), expression(mu))
segments(x0 = 2, y0 = mean(log10(doubling$r)[doubling$kingdom == "Bacteria"]), x1 = 2.39, y1 = mean(log10(doubling$r)[doubling$kingdom == "Bacteria"]), lwd = 2, lty = 2, col = "black")
text(x = 2.45, y = mean(log10(doubling$r)[doubling$kingdom == "Bacteria"]), expression(mu))
legend(x = 0.3, y = 0.8, paste("n =", length(doubling$species), sep = " "), cex = 0.9, bty = 'n')
legend(x = 0.3, y = 0.6, paste("p < 0.001"), cex = 0.9, bty = 'n')
legend(x = 0.3, y = 0.4, expression(R^2==0.019), cex = 0.9, bty = 'n')
#mean(doubling$doubling[doubling$kingdom == "Archaea"])
#mean(doubling$doubling[doubling$kingdom == "Bacteria"])
stripchart(log10(r) ~ kingdom, vertical = TRUE, data = doubling, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)
mtext("A", side = 2, at = 1, line = 2.8, las = 1, font = 2)
## Basic boxplot
vioplot(d1_up ~ kingdom, data = d1_up, border = NA, method = "jitter", side = "right", ylab = expression(paste("Cell size")~log[10]~(mu*m)), xlab = "Kingdom", col = "white", las = 1)
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
vioplot(d1_up ~ kingdom, data = d1_up, border = NA, method = "jitter", side = "right", ylab = "", xlab = "", col = c(alpha("purple", 0.2), alpha("orange", 0.2)), las = 1)
segments(x0 = 1, y0 = mean(d1_up$d1_up[d1_up$kingdom == "Archaea"]), x1 = 1.311, y1 = mean(d1_up$d1_up[d1_up$kingdom == "Archaea"]), lwd = 2, lty = 2, col = "black")
text(x = 1.4, y = mean(d1_up$d1_up[d1_up$kingdom == "Archaea"]), expression(mu))
segments(x0 = 2, y0 = mean(d1_up$d1_up[d1_up$kingdom == "Bacteria"]), x1 = 2.4, y1 = mean(d1_up$d1_up[d1_up$kingdom == "Bacteria"]), lwd = 2, lty = 2, col = "black")
text(x = 2.45, y = mean(d1_up$d1_up[d1_up$kingdom == "Bacteria"]), expression(mu))
text(x = 0.65, y = 5.3, paste("n =", length(d1_up$species), sep = " "), cex = 0.9)
text(x = 0.65, y = 4.97, paste('p < 0.001'), cex = 0.9)
text(x = 0.68, y = 4.68, expression(R^2==0.043), cex = 0.9)
#mean(d1_up$d1_up[d1_up$kingdom == "Archaea"])
#mean(d1_up$d1_up[d1_up$kingdom == "Bacteria"])
stripchart(d1_up ~ kingdom, vertical = TRUE, data = d1_up, method = "jitter", add = TRUE, pch = 20, col = c(alpha("purple", 0.3), alpha("orange", 0.5)), cex = 1.3)
mtext("B", side = 2, at = 6, line = 3, las = 1, font = 2)
#dev.off()Model diagnosis
size <- lm(d1_up ~ kingdom, data = d1_up)
summary(size)>
> Call:
> lm(formula = d1_up ~ kingdom, data = d1_up)
>
> Residuals:
> Min 1Q Median 3Q Max
> -1.1004 -0.2783 -0.0783 0.1217 4.6217
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 1.25044 0.04059 30.803 <2e-16 ***
> kingdomBacteria -0.37217 0.04339 -8.577 <2e-16 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 0.5727 on 1594 degrees of freedom
> Multiple R-squared: 0.04412, Adjusted R-squared: 0.04352
> F-statistic: 73.57 on 1 and 1594 DF, p-value: < 2.2e-16layout(matrix(c(1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
par(las = 1)
plot(size)
growth <- lm(doubling_h ~ kingdom, data = doubling)
summary(growth)>
> Call:
> lm(formula = doubling_h ~ kingdom, data = doubling)
>
> Residuals:
> Min 1Q Median 3Q Max
> -13.64 -10.87 -8.89 -1.36 421.12
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 13.867 2.071 6.695 3.75e-11 ***
> kingdomBacteria -1.981 2.380 -0.832 0.405
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Residual standard error: 31.07 on 926 degrees of freedom
> Multiple R-squared: 0.0007477, Adjusted R-squared: -0.0003314
> F-statistic: 0.6929 on 1 and 926 DF, p-value: 0.4054layout(matrix(c(1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4), nrow = 4, ncol = 4, byrow = TRUE))
par(las = 1)
plot(growth)
Figure S1
## Relationship between diameter and volume
d1_up <- aggregate(spp.d1_up$d1_up, by = list(spp.d1_up$species), mean, na.action = na.rm)
cell.vol <- aggregate(vol$volume, by = list(vol$species), mean)
names(d1_up) <- c("species", "d1_up")
names(cell.vol) <- c("species", "volume")
d1_up$species <- gsub("[[:punct:]]", "", d1_up$species)
d1_up$species <- gsub(" ", "_", d1_up$species)
cell.vol$species <- gsub(" ", "_", cell.vol$species)
di.vol <- merge(cell.vol, d1_up, by = "species")
p.cor <- cor.test(log10(di.vol$volume), di.vol$d1_up, method = "pearson")
##png("figure3.png", height = 7, width = 7, units = "in", res = 360)
plot(log10(volume) ~ d1_up, data = di.vol, pch = 16, las = 1, ylab = expression(paste("Cell volume ", log[10], sep = " ")(mu*m^3)), xlab = expression(paste("Cell diameter ", log[10], sep = " ")(mu*m)), mgp = c(2.6, 1, 0), type = "n")
grid(nx = NULL, ny = NULL, col = alpha("lightgray", 0.5), lwd = 1, lty = 2)
par(new = TRUE)
plot(log10(volume) ~ d1_up, data = di.vol, pch = 16, las = 1, ylab = "", xlab = "", mgp = c(2.6, 1, 0), type = "p", axes = FALSE)
abline(lm(log10(di.vol$volume) ~ di.vol$d1_up), lwd = 2)
legend(x = 0.1, y = 2.3, legend = expression("R"^2==~0.877), bty = "n")
legend(x = 0.1, y = 2.06, legend = expression(p==~0.001), bty = "n")
##dev.off()